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SUMMARY 


TWO-BOUNDARY GRID GENERATION FOR THE SOLUTION OF THE 
THREE-DIMENSIONAL COMPRESSIBLE 
NAVIER-STOKES EQUATIONS 


Robert Edward Smith 


A grid generation technique called the "two-boundary technique" is 
developed and applied for the solution of the three-dimensional com- 
pressible Navier-Stokes equations describing laminar flow. The Navier- 
Stokes equations are presented relative to a xyz cartesian coordinate 
system and are transformed to a Cn? computational coordinate system. 
The grid generation technique provides the Jacobian matrix describing 
the transformation. 

The "two-boundary technique" is based on algebraically defining 
two distinct boundaries of a flow domain and joining these boundaries 
with a "connecting function" which is proposed to be linear or cubic 
polynomials. The algebraic boundary representation can be analytical 
functions or numerical interpolation functions. Control of the distri- 
bution of the grid in the physical domain is achieved by embedding "con- 
trol functions" which redistribute the uniform grid of the computational 
domain and concentrate or disperse the grid in the physical domain. The 
computer program to solve the Navier-Stokes equations is based on a 



MacCormack time-split technique and is specifically designed for the 
vector architecture and virtual memory of the CYBER 203 computer. The 
program "Navier-Stokes solver" is written in the SL/1 language which 
allows 32-bit word arithmetic operations and storage. The program can 
run with 5 x 10^ grid points using only primary memory, and the compu- 
tational speed is 4 X 10"^ seconds per grid point per time step. 

Using the "two-boundary technique," grids are developed for two 
distinctly different flow field problems, and compressible supersonic 
laminar flow solutions are obtained using the Navier-Stokes solver. 
Grids and solutions are obtained for a family of three-dimensional 
corners at Mach number 3.64 and Reynolds numbers 2.92 x lO^m and 
3.9 X 10^/m. Also, grids are derived for spike-nosed bodies, and solu- 
tions are obtained at Mach number 3 and Reynolds number 7.87 x 10^/m. 

Coupled with the Navier-Stokes solver, the "two-boundary technique" 
is demonstrated to be viable for grid generation associated with com- 
puting supersonic laminar flow. The technique is easy to apply and is 
applicable to a wide class of geometries. The "two-boundary technique" 
can serve as the foundation for generating grids with highly complex 
boundaries and yield grid point distributions that can capture rapidly 
changing variables in a flow field. 
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1. INTRODUCTION 

In recent years, the availability of large scale scientific com- 
puter systems has resulted in rapid progress in the field of Computa- 
tional Fluid Dynamics. There is now the capability to calculate many 
complex unsteady two-dimensional and steady three-dimensional flows. 
MacCormack and Lomax [1]* summarize the "state of the art" for the 
computation of compressible viscous fluid flow. For a heat conducting 
compressible fluid acting near body surfaces with large separation 
regions or inviscid-viscid interactions, the numerical solution of the 
Navier-Stokes equations is the preferred approach [1]. An emerging 
problem, however, is the generation of grid systems on which solutions 
can be obtained when there are complex boundary geometries. This prob- 
lem is compounded in three dimensions. This study addresses the solu- 
tion of the three-dimensional compressible Navier-Stokes equations, 
the generation of grids, and the solution algorithm-computer relation- 
ship. The emphasis is placed on grid generation. 

An algebraic grid generation technique applicable to the Navier- 
Stokes equations is developed, and a three-dimensional Navier-Stokes 
solver (compressible laminar flow) based on a proven numerical tech- 
nique (MacCormack time-split algorithm [1-4]) is developed for the CDC 
CYBER 203 vector computer [5]. Also, flow visualization techniques have 
been developed in conjunction with this research but will not be dis- 
cussed in detail. In order to evaluate the overall system for computing 
viscous compressible flow, and in particular the grid generation 


*The numbers in brackets indicate references. 
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technique, grids are determined for a family of three-dimensional 
corners and two spike-nosed bodies. • 

The grid generation technique is called the "two-boundary tech- 
nique." It is applicable in two and three dimensions and is a method- 
ology for direct computation of the physical grid as a function of a 
uniform rectangular computational grid. The Jacobian matrix of the 
transformation can be obtained by direct analytic differentiation. This 
is in contrast to the indirect approach where an elliptic partial dif- 
ferential equation system is solved for the coordinates of the physical 
grid relative to the computational grid, and in which the Jacobian 
matrix must be obtained by numerical differentiation. The indirect 
approach is popularly known as the "TTM method" [1, 6-10]. In the 
"two-boundary technique," two separate non-intersecting boundaries are 
defined by means of algebraic functions or numerical interpolation 
functions. These functions have as independent variables, coordinates 
which are normalized to unity. Another function with an independent 
variable defined on the unit interval connects the boundaries. 

The "two-boundary technique" is based upon concepts found in the 
theory of surface definition [11,12]. Gordon and Hall [13] postulate 
the essentials of the technique and emphasize finite element grids. 

Also, Eiseman [14-16] uses a form of the technique in generating grids 
for multi connected two-dimensional domains. In this investigation the 
"two-boundary technique" is developed and is analyzed for finite differ- 
ence solutions for fluid flow applications. Low order polynomials 
(linear and cubic) are used for connecting functions. For the cubic 
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connecting function, orthogonality can be enforced at the boundaries 
through knowledge of the normal derivatives there. Control of the grid 
(grid spacing in the physical domain) is achieved by the superposition 
onto the independent variables algebraic or transcendental functions 
with desirable characteristics. Splines under tension [17-19] are pro- 
posed for approximate boundary definition. The "two-boundary technique" 
is used to algebraically generate grids for a family of three-dimensional 
corners and to generate a combined algebraic-numeric grid for spike- 
nosed bodies. The derivatives composing the Jacobian matrix for the 
three-dimensional corners and spike-nosed bodies are presented for 
obtaining numerical solutions of the Navier Stokes equations. 

The CDC CYBER 203 is a large scale computer with vector processing 
architecture and virtual memory. Generally efficiency using a vector 
computer increases with increasing vector length, however, considerable 
attention must be given to the algorithm-machine architecture relation 
and balancing the vector length with practical limits of primary memory. 

A MacCormack time-split solution algorithm is programmed for the 
CYBER 203 computer and is called the "Navier-Stokes solver." The 
MacCormack technique is used because of its robustness and adaptability 
to vector processing. Another primary consideration when developing a 
"Navier-Stokes solver" on a large complex computer is the capability to 
solve a wide class of problems with a minimum of programming changes. 

This has been accomplished by programming the complete transformed 
equations of motion and storing all nine elements of the Jacobian 
matrix of the transformation at each grid point (transformation data). 
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Supplying the transformation data from a grid generation technique and 
programming the boundary conditions 'for a given problem (separate sub- 
routine) allows the program be applied to virtually any laminar fluid 
flow problem. Since the split MacCormack technique is used, two- 
dimensional solutions can be obtained without unnecessary computations. 
The operator for the third dimension is bypassed. A final important 
point relative to the Navier-Stokes solver is that the MacCormack tech- 
nique is written in the SL/1 language [20] and uses the 32-bit arithmetic 
option of the CYBER 203. By using 32-bit words, twice the in-core stor- 
age is available and approximately twice the computational speed is 
achieved compared to the use of normal 64-bit words. There are approxi- 
mately two million words of primary memory and the computational speed 

_5 

is 4 X 10 seconds per grid point per time step for the 32-bit word 
length. For the explicit technique, no significant degeneration in 
accuracy is observed using the smaller word size. The Navier-Stokes 
solver is independent of the grid generation technique, and the trans- 
formation data from any technique can be used by the code. 

Using the "two-boundary technique" grids are developed for two 
distinctly different flow field problems, and compressible supersonic 
laminar flow solutions are obtained using the computer program based on 
the MacCormack technique. A set of algebraic grid generation equations 
are developed using the "two-boundary technique" for a family of three- 
dimensional corners consisting of wedge-cylinder, plate-cylinder, 
approximate wedge-plate, and approximate rectangular corners. It is 
also shown that exact grids for planar intersecting corners can be 
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derived with the "two-boundary technique." Corner flow solutions are 
obtained on a 20 x 36 x 36 grid and a 12 x 64 x 64 grid. The solutions 
obtained on the 12 x 64 x 64 grid are compared with physical experiments 
and other numerical experiments. The Mach number used is 3.64 and the 
Reynolds number is 2.92 x 10^/m and 3.9 x 10^/m. 

Also, algebraic grids are derived using the "two-boundary technique" 
for spike-nosed bodies. In particular, grids for a one-half inch spike- 
nosed body and a one and one-half inch spike-nosed body are obtained. 
Supersonic flow solutions at Mach number 3 and Reynolds number 
7.87 X 10®/m are obtained about these configurations. Unlike the flows 
about the three-dimensional corners, the flow about the spike-nosed 
bodies is unsteady. The amplitude of the oscillations about the one-half 
inch nose body is quite small, however, the one and one-half inch spike- 
nosed body flow field oscillates with a large amplitude. The high 
amplitude solutions are compared with physical experiments. The flow 
fields are two-dimensional axisymmetric, but are solved with a three- 
dimensional Navier-Stokes solver resulting in considerable savings of 
development time for a specialized axisymmetric code. 

For flow visualization, a relatively novel approach has been 
developed where a color spectrum is used to display a scalar variable 
such as density, Mach number, etc., on a two-dimensional slice of a flow 
field. Sequences of pictures can show the history of a developing flow 
or a scan of the flow field in a three-dimensional domain. The Diccomed 
Digital Display/Film Writer system which is normally used for environ- 
mental image processing is used for the flow visualization. 
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In summary, the main objectives of this study are the development 
of an algebraic grid generation procedure, the development of software 
to solve the compressible three-dimensional Navier-Stokes equations on 
a vector computer using the results of the grid generation technique, 
and the application of the grid generation technique and software to 
solve specific supersonic flow problems. The organization is as 
follows. In Chapter 2 the three dimensional compressible Navier-Stokes 
equations are presented relative to a Cartesian coordinate system and 
are transformed to a uniform grid computational coordinate system. 

This introduces the information that must be determined by the grid 
generation technique. The "two-boundary technique" is developed and 
applied to generate grids and Jacobian derivatives for a family of 
three-dimensional corners, spike-nosed bodies, and an airfoil configura- 
tion. In Chapter 3, the MacCormack technique is presented, and its 
compatibility with the CYBER 203 is described. In Chapter 4, supersonic 
flow solutions about three-dimensional corners and spike-nosed bodies 
obtained with the "two-boundary technique" and Navier-Stokes solver 
are described. 

2. ANALYSIS 

This chapter develops the equations of motion and the "two- 
boundary technique" for grid generation. Grids and boundary conditions 
are developed for a family of three-dimensional corners and for spike- 
nosed bodies. Also, grids are developed for airfoil boundaries using 
splines under tension. 
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2.1 Navier-Stokes Equations of Motion 
The governing equations which describe the motion of a viscous 
compressible heat conducting fluid are the continuity equation, momen- 
tum equations, and energy equation. These equations are derived from 
the concept of continuum mechanics. The continuum concept and deriva- 
tion of the Navier-Stokes equations of motion are found in several 
references, of which Schlichting [21] is the most notable. 

Expressed in symboic form the Navier-Stokes equations of motion 

are: 


Continuity: 

If + V • (pID = 0. 

(2.1a) 

Momentum: 

V • (pGQ - t) = 0, 

(2.1b) 

Energy : 

+ V • (peu + q - u • t) - 0. 

(2.1c) 


The stress tensor, dissipation function, and heat conduction for a 
rectangular cartesian coordinate system are: 

= stress tensor 
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function of temperature and is 
semi empirical equation: 


is set equal to zero. This is 
gas where the molecules has no 
internal degree of freedom. For a polyatomic gas the bulk viscosity is 
not always zero and can be the same order of magnitude as the molecular 
viscosity in sound propagation and shock structure. A detailed discus- 
sion of bulk viscosity is given by Vincenti and Kruger [22]. 

At this point there are five coupled partial differential equations 
and one algebraic equation with eight unknowns: p, u, v, w, P, e, 

T, and y. In order to have a complete system, there must be two addi- 
tional equations relating the unknowns. The equation of state is 

P = P (P,T), 

and for a perfect gas P = pRT and e = C^T where Cy is the specific 
heat at constant volumn, and R is the gas constant. For compatible 
boundary conditions this system of equations is solvable. 


The viscosity coefficient y is a 
adequately approximated by Sutherland's 


with 


y 

Pv. 


T + S 

^ 

T + S_ 



S = 198.6°R for air. 
0 


The bulk viscosity coefficient y^ 
a reasonable assumption for a monatomic 
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2.2 Transformed Equations of Motion 
The equations of motion in Section 2.1 are expressed in terms of a 
Cartesian coordinate system. If an object is defined in this coordinate 
system and a flow is to take place about the object, it is desirable to 
perform the computation in a coordinate system which conforms to the 
boundaries of the object. There are two primary reasons for wanting the 
coordinate system to be boundary-fitted. Boundary-fitted coordinates 
afford the ability to apply boundary conditions exactly avoiding inter- 
polation error, and they minimize the logic that is necessary to apply 
boundary conditions. The penalty for these advantages is added com- 
plexity of the equations of motion. Another consideration is that when 
the domain of the flow field is discretized, it is desirable to have 
grid points concentrated in certain regions where high rates of change 
are likely to occur. For instance, in the boundary layer region more 
grid points are necessary to resolve the rapid change in the state 
variables. If the cartesian coordinate system where the object is 
defined is called the physical domain, the coordinate system relative 
to the boundaries of the object is called the computational domain. The 
relationship between the physical domain and the computational domain is 
a unique single-valued transformation with continuous derivatives such 
that if the coordinates in the computational domain are 5, n> C: 
then 


? = ?(x,y,z), n = n(x,y,z), and ? = c(x,y,z). 
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Conversely, 


X = x(5,n.c), y = y(C.n.c), and z = z{5,n,c) 

where x, y, and z are coordinates in the physical domain. Since 
the equations of motion in terms of the Cartesian coordinate system of 
the physical domain are advantageously solved in terms of the coordinate 
system of the computational domain, the equations must be transformed. 
This is accomplished by expressing the derivatives of the state variables 
with respect to the components of the physical domain in terms of 

the CnC components of the computational domain as follows: 
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Notice that u, v, and w are the velocities along the x, y, 
and z axes in the physical domain. 
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is called the Jacobian matrix of the transformation whose elements are 
the nine derivatives specifying the rate of change of the computational 
coordinates with respect to the physical coordinates. The equations of 
motion (Eq. (2.1)) are in conservation form [3] and can be expressed as 


8U . 3F . 3G . 3H _ n 

H ^ ST 57 3? ' “ 


( 2 . 2 ) 


where 


U 


— - 


pu 

p 

pu 


PUU - T 

XX 

pv 

F = 

PUV - T 

xy 

pw 

pe 


Puw - T 

X2 

f 



peu + q^ - 


— 

pv 


* 

pw 

PUV - 


puw - 

pw - 

, H = 

pVW - 

PVW - 


pww - 

pev + % - 


pew + - <(>2 _ 
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Equation (2.2) written with respect to the transformed variables of the 
computational domain becomes: 



K E. 

9? 9x 


+ 


9F 

M + 

3F 

3? 

+ ^ 

E + 

3n 

9x ^ 

3? 

9x 

^ 3? 

3y 


9H 

E + 

9H 

E + 

9H 

3C _ 

3 ? 

9z ^ 

3n 

92 ^ 

3? 

9z 



9? 9y 


(2.3a) 


This can be written compactly as: 



(2.3b) 


where 


"l 


K 


^1 


X 


^1 


F 

^2 

= 

n 

9 

^2 


y 

9 

^2 

= 

G 

_^3_ 


c 


_^3_ 


z 


/3_ 


H 


Given the Jacobian matrix at each grid point and initial and boundary 
conditions, the transformed governing equations of motion are in a form 
to be numerically solved. It is noted at this point that the equations 
of motion are in weak conservation form relative to the Jacobian 


deri vatives. 
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2.3 Definition of a Computational Domain and Transformation Data 
In the previous section the Navier-Stokes equations are transformed 
from a Cartesian coordinate system to a computational domain. In so 
doing nine additional unknowns, which are the elements of the Jacobian 
matrix, are added to the problem. When the finite difference technique 
described herein is applied to Equation (2.3), the Jacobian matrix 
must be known at each grid point. The objective of a grid generation 
technique is to provide the Jacobian matrix which is henceforth called 
transformation data. The computational domain is defined in this sec- 
tion along with the formulas necessary for computing the transformation 
data based on known functional relations between the computational domain 
and the physical domain. The next three sections concentrate on deter- 
mining functional relations between the computational and physical 
domains. The computational domain is defined to be a rectangular 
parallelepiped and a uniform grid is superimposed onto the domain 
(Fig. 1) such that: 

0 < ? < 1 , 

0 £ r) £ 1 , 

0 £ ? £ 1 , 

(2.4) 

A? = constant-j , 

An = constant2» 


Ac = constant^. 
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A functional relation between the computational domain and the physical 
domain can be expressed as 

X = x(5,n,?), y = y(C,n.?). and z = z(^,n,c). (2.5) 

Further, these functions must map boundaries in the computational 
domain onto boundaries in the physical domain such that 

Xg = x(5B.nB.CB)» Vb = 

where Xb» yB and Zb define the boundaries of the physical domain 
and Cb> ^B boundaries of the computational 

domain. The transformation data is composed of the rates of change of 
the computational coordinates with respect to the physical coordinates. 
If the inverse functional relations 

? = ^(x,y,z), n = n(x,y,z), and ? = c(x,y,z) (2.6) 

are known, the transformation data can be directly found by differ- 
entiation. It is not necessary, however, to know the inverse functional 
relations to determine the transformation data. The Jacobian matrix 
can be evaluated by differentiating the functional relation (Eq. (2.5)). 
That is 


3X 


3x 

3? 

3n 

3C 

9y. 


3^ 

3C 

3n 

3C 

9Z 

31 

3z 

3? 

3n 

3? 
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and then 


-1 


j _ Transposed of Cofactor (J~ ) 


J 


-1 


5 


where |J is the Jacobian determinate and J"^ is the inverse 
Jacobian matrix. 


3X 

3x 

3x 

9^ 

9n 

9? 

h. 



95 

9n 

9C 

3z 

3z 

3z 

95 

9n 

9? 


and 


_ 9x 

(1^ 

3z 

_ ^ 

— ) 
9n^ 

3x 

95 

^9n 

9? 

3? 

" 9n 

. 3x 

^9C 

3Z 

. 

92) 

95^ 


9C 

9n 

9n 



3Z. ^ 

3? ■ 3? 3^^ 


95 

95 

95 


3x 


3X 

3x 

9y 

3Z 


95 

3n 

9? 

9n 

3n 

9n 

= J = 



9Z 

3X 

9y 

3Z 


95 

9n 

9C 

9C 

9? 

9C 


3z 


3z 

3x 

9y 

3Z 


9? 

3n 

9? 
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^3n 3C " 3C 3n^ 


/ 3x ^ ^ 3Z x 

3n 3? 3? 3i1 


/ 3x 3y 3x 3 Y n 
^3n 3? ■ H 3n^ 


J = 


U'^l 


^35 3^ ■ 3? 3?^ 

(^ il. 3^ ^ 

^3? 3n " 3n 3C^ 


/ 3x 3z 3x 3Z n 

^35 3C “ 3? 35^ 

/ 3x 3z 3x 3z ^ 

^3? 3n “ 3n 3C^ 


/ 3x 3y 3x 3y x 

^35 3C ~ 3? 35^ 

^3^ 3n ■ 3n 35^ 


(2.7) 

provided |J”^| ?* 0. 

The transformation data can be pre-evaluated and stored or it can 
be computed as needed. The trade off is the additional computation 
cost versus the storage cost. For the Navier-Stokes solver discussed 
in this study, the transformation data is precomputed and stored for 
later use. 


2.4 Two-Boundary Grid Generation 
A computational domain is postulated by Equation (2.4). It is a 
rectangular parallelpiped in three dimensions and a square in two 
dimensions. The physical domain is a subdomain of a Cartesian coordi- 
nate system. A transformation between the physical domain and the com- 
putational domain is a mathematical relationship mapping one domain 
onto the other. Similarly a grid in one domain is mapped onto a grid 
in the other domain. When the transformation maps boundaries in the 
physical domain onto boundaries in the computational domain the term 
"boundary-fitted coordinate system" is used to describe the 


transformation. 
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An indirect (differential) approach for finding the relationship 
between the computational and physical grids described by Thompson 
et al. [6-10] has been highly successful. In this approach the 
elliptic system of partial differential equations which must be satis- 
fied by the mapping between the two domains is numerically solved by an 
iterative technique such as Successive-Over-Relaxation (SOR). The 
numerical solution is the grid in the physical domain corresponding to 
the grid in the computational domain. The transformation data is 
obtained by numerical differentiation, and a grid change requires a 
new solution of the elliptic system. 

A direct (algebraic) approach, where an explicit functional rela- 
tionship between the computational domain and the physical domain is 
known, has the advantages that changes to the grid are direct, rapidly 
obtained, and transformation data is analytically available. 

A direct algebraic approach called the "two-boundary technique" 
is described in the present paper. The technique has a wide variety 
of applications in both two and three dimensions. A preliminary 
description of the technique is presented in [23]. Symbolically, the 
relation between the computational domain and physical domain can be 
written as 


X = X(?,n,?), (2.8a) 

y = Y(C,n.?), (2.8b) 


^ 2(^,n» c) , 


(2.8c) 
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0 < ? < 1 , 

0 £ n £ 1 . 

0 £ c £ 1 . 

Equation (2.8) is equivalent to Equation (2.5). For a boundary-fitted 
relationship between the two domains, boundaries in the computational 
domain should map onto boundaries in the physical domain as shown in 
Figure 2. For instance, for the boundaries q = 0 and n = 1 in the 
computation domain. Equation (2.8) becomes 

= X(^,0,c) = x^(^,c), 
y^ = Y(C,0,?) = Y^(?,?), 

z-j = 2(^,0,^) = Zi(^,?), 

^2 ~ X(^,l , 5 ) - i 

y2 “ Y(^,l,?) = Y2 (C,c), 

Z2 = 2 ( 5 , l,c) = 22(5,5). 

Here, X-|(5,5)» X2(5»5), etc. are boundaries in the physical domain 

and, as such, are functions defined only at the boundaries. An approp- 
riate explicit expression for Equation (2.8) would separate one 
variable (q) to be independently varied with parameters dependent on 
position and derivatives on the boundaries. Since the boundaries are 
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themselves functions and can be determined independently. Equation (2.8) 
can be rewritten as 


X = X(c,n,c) = x(x^(?,?). X2(?.c), n), 

(2.9a) 

y = Y(^.n,?) = Y(Y^(e,c). Y2 (c, 0. n). 

(2.9b) 

z = Z(^,n,c) = Z(z^(^,i;), n). 

(2.9c) 

The explicit forms of Equation (2.9) proposed herein 
metric linear and cubic polynomials. 

Linear 

are simple para- 

X = X2(?,?)n + X^(c,?)(l - n), 

(2.10a) 

y = Y2(C,c)n + Y^(c,?)(l - n), 

(2.10b) 

z = Z 2 (?,c)n + Z^(?,c)(l - n). 
0 £ n < 1 » 

(2.10c) 


Cubic 

, dX, 

X = x^(?.c)fi(n) + X2(c,c)f2(n) + — ^ (5.df3(n) 

dn 

dXp 

dn ^ 


(2.11a) 
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UI-. 

y = Y^(?.c)fi(n) + Y2(^,c)f2(n) + (C.?)f3(n) 

dY2 

+ — (?,?)f4(n). 

dn ^ 


(2.nb) 


U^-| 

z = Z,(c,?)fi(n) + Z2(C,?)f2(n) + — (?,?)f3(n) 

dn 

dZ. 

+ -^(?.c)f4(n), (2.11c) 

dn 


where: 

f-| (n) = 2n^ - 3n^ + 1 , 

f2(^) “2n^ + 3n^, 

f3(n) = n^ - 2n^ + n> 

■p4(n) = n^ - n^. 

0 ^ n ^ 1 • 

A function such as Equation (2.10) or Equation (2.11) is topologically 
referred to as a homotopy [24]. Blending-function [11] is another 
name that has been given to such equations for problems in surface 
design. Herein, because of the context in which they are used, they 
are defined as "connecting functions." 
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Applying a cubic connecting function implies that the physical 

grid can be forced to be orthogonal' at the boundaries since the deriva- 
dX, dXp 

tives etc. can be computed from the cross product 

dX, dXi dY, dY, 

of the tangential derivatives 
etc. That is, 


^ dY, ^ dZ, 

^ (?.C) t + — (C.?)J + — (?,c)k = 


dn 


dn 


dn 


-»• 

J 


t 




d? de 


dC 




>P dY„ dZ 

— (?.?) — (?.C) —(?,?) 


d? 


d? 


d? 


£ = 1,2 


where i, j, and k are unit vectors and K is the magnitude of the 
normal vector. Applying this procedure will force the grid to be 
orthogonal at the boundaries but not necessarily anywhere else. For the 
linear connecting function, the physical grid will seldom be orthogonal. 

Given the connecting function and parametric boundary functions, 
a uniform computational grid can be mapped onto the physical domain 
forming a physical grid. Concentration of grid points in the n direc- 
tion is accomplished by choosing a function ri = fi(ri) such that 
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0<n<l» 0<fi<l, and '^>0 (Fig. 3). For example, contracting 
the physical grid towards one boundary or the other can be accomplished 
by 


n = ; 0<n<l. (2.12) 

e*^ - 1 

where k is a free parameter whose magnitude dictates the degree of 
contraction. Embedding this exponential function in the linear con- 
necting function, Equation (2.10) becomes: 

X = X 2 (C,c)n + X^(?,?)(l - n), 

y = Y2(5,?)n + Yi(5,c)( 1 - n), 

z = Z 2 (?,?)n + Zi(C,c)(l - n), 

0 £ n j< 1 . 

Once the connecting function has been chosen, the remaining prob- 
lem is the determination of the boundary functions which are independ- 
ent of n. For the "two-boundary technique" the approach is to choose 
parametric variables i and t associated with the boundaries such 
that 

x-|(5,?) ->■ X^(s,t), 


yi(?jC) Y-j(s,t), 
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z^(?,C) -»■ Z^(s,t), 


X2(5.?) X2(i,t) 

y2(5.?) -*■ V2(s,t), 
Z2(C,t) -> Z2(s,t) 


0 <£< 1 , S. <^<S , 

- ^ ’ min - - max’ 


0<C<1, t. <t<t 
- ^ - ’ mm - - max 


The choice of parametric variables can vary from problem to problem. 
A relationship between (C.c) and (s,t) is 


s = C(s* - s . ) + s . , 
max min' min’ 


t = c(t “ t • ) + t . . 
max min^ min 


This is a linear relation which maps the unit interval onto the 
parametric variables. Control of the physical grid at the boundaries 
is accomplished in the same manner as for the connecting function. 




c = c(c). §> 0, 



Since the connecting function is dependent on the boundary position, 
control of the entire grid is accomplished. 

2.4.1 Approximate Boundary-Fitted Coordinate Systems 
Using Tension Spline Functions 

It is often the case that boundaries in a physical domain are 
described by discrete sets of points. The boundaries may be open or 
closed (Fig. 2). An approximate boundary-fitted coordinate system can 
be obtained using the "two-boundary technique" and a tension spline 
function interpolation to the discrete data defining the boundaries. 
Tension splines [17-19] are chosen because standard cubic splines [25] 
and other higher ordered interpolation techniques often result in 
wiggles in the approximation. Wiggles on a boundary using the "two- 
boundary technique" propagate into the interior grid. The tension 
parameter embedded in the tension spline interpolation allows control 
of the "curvedness" of the approximation. A very large magnitude of 
the tension parameter corresponds to a linear interpolation whereas a 
very small value corresponds to cubic splines. Tension splines can 
be applied in two and three dimensions. A two-dimensional example is 
presented. 

Using the tension spline technique, a point set on boundary one 
is defined by {x^. ,y^. and on boundary two by {Xj ,yj}j~^|’. 
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Approximate arc length is used as a parametric independent variable. The 
approximate arc length is: 

V 2^/2 

S, = [(X,^., - X,) -y,)] 

2 2 

i = 1 . . . n 
j 1 « • • m 

s, = 0 

0 1 s, < s„ 

0 < s. < s . 

— j — m 

From the computational coordinate system the unit interval (0 £ C £ 1) 
must be mapped onto each boundary; that is: 

s = s(?), 

0 £ 5 £ 1 . 

This is accomplished by letting 

s = ^s^ on boundary one and 
s = on boundary two. 
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The tension spline function is a piecewise continuous set of trans- 
cendental functions where x and y between the I and + 1 points 
are defined by 

sinh[a(sj , - s)] 

X = g"(s ) 

a sinh[a(s^^^ - s^)] 

sinh[a(s - s )] 

“2 

aSinh[a(s^^^ - s^)] 




y = h"(s^) 


sinh[a(s^^^ - s)] 
a^sinhC (Sj^, - Sj)] 




sinh[a(s - s^)] 
a^sinh[a(s^^^ - s^)] 


( 2 . 13 ) 
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(2.14) 


s - s(C) ^’S^ax’ 

5, = i on boundary one, 
£ = j on boundary two, 
a = tension parameter. 


The unknowns in these equations are g"(s£) and h"(s^) which 

£=N 

are second derivatives at the data points {x£,yj ^}^_1 where N = n for 
boundary one, N = m for boundary two, and are obtained through enforce- 
ment of the continuity of the first derivatives at the data points and 
the specification of two end conditions. A tridiagonal system of 
linear equations results for each set of unknowns. The solution of the 
tridiagonal systems yield g"(s£) and h'Hs^). 

The cubic connecting function (Eq. 2.11), and the exponential 

function (Eq. 2.12) provide the relationship between the computational 

dYo 

domain and the physical domain. The derivatives and 


and 


are: 
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= K ^ 

dn 

ds 



dn 

ds 


By defining a grid with constants and An in the computational 
domain a corresponding grid is explicitly defined in the physical 
domain. 

An example of a grid about a Kcirm^n-Trefftz airfoil is presented 
using the spline under tension approximation to the boundaries and a 
cubic connecting function. Table 1 contains the data describing the 
airfoil boundary and outer boundary. Figure 4 shows the approximation 
to the airfoil boundary and Figure 5 shows the grid. A tension parame- 
ter value of 2 is used. Transformation data have not been computed 
for this example. 

2.4.2 Transformation for a Wedge-Cylinder Corner 

An application of the two-boundary technique using analytical sur- 
face functions is a family of three-dimensional corner geometries which 
occur in many aerodynamic situations (Fig. 6). Supersonic flow about 
these geometries is characterized by strong visid-inviscid interactions 
which are adequately analyzed only through the numerical solution of 
the Navier-Stokes equations. When solving this system of equations 
with a finite difference technique, a grid must be designed to capture 
the interactions and allow for accurate application of the boundary 
conditions. 



Table 1. Data description for an airfoil grid 


Inside boundary 

X ft. 

y ft. 

.49950 

-.000031 

.49860 

-.001400 

.49600 

-.002760 

.48620 

-.005550 

.47060 

-.008510 

.39010 

-.028590 

.26960 

-.029970 

.12270 

-.040790 

-.03480 

-.048450 

-.18750 

-.050520 

-.32110 

-.045590 

-.42390 

-.0377390 

-.48650 

-.016530 

-.50270 

.003820 

-.47110 

.026640 

-.39690 

.048640 

-.28790 

.066100 

-.15380 

.075290 

-.00560 

.074530 

.14400 i 

.064350 

.28230 

.04750 

. 39580 

.028260 

.47200 

.011300 

.48710 

.006810 

.49510 

.003040 

.49880 

.001430 

.49950 

-.000031 


Outside boundary 

X ft. 

y ft. 

0.0 

3.0 

2.12 

2.12 

3.0 

0.0 

-2.12 

-2.12 

0.0 

-3.0 

-2.12 

-2.12 

-3.0 

0.0 

-2.12 

2.12 

0.0 

3.0 






Fig. 5 Grid for Karman-Trefftz airfoil obtained with the "two-boundary technique". 


OJ 

00 
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The "two-boundary technique" is applied to the wedge-cylinder cor- 
ner with the aid of Figure 7. The other corner geometries are derived 
from the wedge-cylinder definition. The physical domain is the region 
enclosed by the circular cylinder with radius rQ, an outer surface 
defined by the wedge angle and a second cylinder radius, and two planes. 
The left plane (wedge surface) is oriented at an angle (p (wedge angle) 
with the longitudinal axis of the cylinder but parallel to the vertical 
axis. The right plane (symmetry plane) is oriented with angle 02 
relative to the vertical axis of the cylinder and includes the longi- 
tudinal axis. The upstream and downstream boundaries are cross sections 
of the region defined by x = Xq and x = Xl and are perpendicular to 
the longitudinal axis. The "two-boundary technique" is applied to this 
geometry by considering the inside cylinder surface as boundary one and 
the outside surface as boundary two. It is desired that n» and ? 
map into the region described above and that 

0 < 5 < 1 , 

0 ^ n £ 1 > 

0 £ ? £ 1 . 

The boundaries are defined by 

rQ,r-j ,Xq,X|^ , 02 >({>» ki ,l<2 
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Boundary one: 

X(?,0,?) = m) = (Xl - Xg)? + Xq, 
Y(C,0,c) = Y.|(5,c) = Tq cos 
Z(?,0,?) = Z^(?,C) = Tq sin 


Boundary two: 

“ X(^) = (Xj^ - Xq)^ + Xqj 

Y(C, 1 ,C) = Y2(C,c) = r^U) cos 
Z(?,l,^) = Z2(^,?) = r^U) sin 

Where: 


= 502 + 0 - 4 ) 01 * 


8, = sin-’ 

‘ \ *^0 > 


03 = sin-^ M^ta riA 


? 2 1/2 

r 2 = [(x(5) tan 


(2.15a) 

(2.15b) 

(2.15c) 


(2.16a) 

(2.16b) 

(2.16c) 
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0 < 5 < 1 , 

0 < C < 1 . 

A linear connecting function is used to generate the internal grid. 
The function is 


X = X2(5,4)n + X^(^,c)(l - n) = X(^) = (x^ - Xq)? + Xq, 

y = Y 2 (?,^)n + Y^(?,00 - n), 

z = Z2(5,?)n + Z^(?,c)(l - n), 


k«ri 

e - 1 



(2.17a) 

(2.17b) 

(2.17c) 

(2.17d) 


0 £ n 1 1 . 

An exponential function is used on both n and ? to concentrate the 
grid in the corner. Figure 8 shows the grid at x = X[_ for cor- 
responding corner surfaces shown in Figure 6. The planar corners are 
closely approximated by letting the radii be very large. 



Rectangular corner 


6 wedge-plate corner 


12.2 wedge-plate corner 



Plate-cylinder corner 6^ wedge-cylinder corner 12.2^ wedge-cylinder corner 

Fig. 8 Grids for wedge-plate and wedge-cylinder corners at x/X|^ = 1. 
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Information needed for the equations of motion (Eq. (2.3)) is the 
transformation data which is obtained from Equation (2.7). The deriva- 
tives in Equation (2.7) are obtained by analytic differentiation of x, 
y, and z (Eq. (2.17)) with respect to C, n> and These deriva- 
tives are; 

^ = X - X 
35 0 



||- = n (5,1,?) + (1 - Ti) If (5.0»?) 

^ = -|Jy(5,o,?) 

If = r| If (5,1,?) + (1 " ^ If (5,0,?) 
|| = n II (5,1 ,?) + (1 - n) H (5,0,?) 
§1 = in zfE 1 cl - 1 cl 

H= n If (5,1,?) + (1 - n) If (5,0,?) 



where 


t (M,C) 
^ (C>1 »c) 

If U,0,?) 
If 

If (?.o,?) 
If (?.l,c) 


_ 99, 

= -r, sinCce, + (1 - C)e,](l - 0 — 

' ^ '95 

9r« _ _ ,_ 

= — cos[502 + (1 - ^) 03 ] - »'2 sin(592 

99, 

- (1 - C) 03 )(l -?) — 

9C 

= -r^ sin[?02 + (1 - C)0 t]C( 02 ' Q])] |f 
= -f 2 sin[c 2 

_ _ 901 

= r, cos[? 05 , + (1 - ?)0-|][l - k ] 

' 95 

9Tj, _ _ _ 

= — sin[c9p + (1 - ^)9,] + rp cos[c02 
9? 

300 

+ (1 - ^goKl - ?) — 

95 

= r^ cos[r02 + (1 - c)0i 3[02 “ 0]3 |f 
= T2 sin[^02 + (1 - c) 03][02 “ 03] |f 
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1 

'7m 

[1 + ] 


9x tan 4> 


^^3 

3C 


1 

TT72 

[1 + [^^] ] 

'^2 


9 /X tan (j)^ 
35 ^ ^2 ’ 


9 /)^ 
35 ^ 


tan 


"■ [r 


9x 

2 tan ()) ^ - X tan 4) 


9r 


35 


hiT^ 



1 X tan 4. II 


The application of this technique in the Navier-Stokes solver is 


found in Chapter 4. 
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2.4.3 Transformation for a Spike-Nosed Body 

The "two-boundary technique" is applied to generate grids about 
spike-nosed bodies (Fig. 9 and Fig. 10). Supersonic flow about these 
bodies is unsteady and separation occurs in the nose-shoulder region. 
Consequently, grids must be concentrated in the nose-shoulder region 
and be adequately spaced to define the shock and the boundary layer on 
top of the shoulder. A linear approximation to the inner boundary and 
a circular arc outer boundary are used in the two-boundary technique. 
Concentration of the grid is accomplished by superimposing an exponen- 
tial function onto the connecting function and a combined exponential 
and parabolic algebraic function is superimposed onto the parametric 
variable along the boundaries. The grid is cast in three-dimensions by 
rotating the two-dimensional description about the axis of symmetry. 

The inside boundary is defined by the set of points 



The parametric variable associated with this boundary is accumulated 
cord length where 


t-| = 0, and 


/V 






,231/2 


/\ 


+ t 


il-1 
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Fig, 10 One and one-half inch spike-nosed body. 
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With the parametric variable defined, two data sets are formed. They 
are: 


I 




A linear approximation to the inside boundary for a spike-nosed body is 
accomplished by linear interpolation of the above data sets. The trans- 
formation data requires the derivatives of the boundary definition with 
respect to the parametric variable and the derivative sets are formed 
and saved for later use. The derivative sets are: 


where 


f 

\ 

M-l 


I 

f 

. dy 

a=M-l 
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t 9 

9 

a 

dt 




“ dt 


1 

a 

i 

|a=l 



a 

a=l 

1— 1 
< X 

-a 


*^a+l 

“X 1 
a- 1 


_ ^a+1 - 

-I 

^a-1 

A. 

dt 

a 


^a+1 

- t 1 

a- 1 

* 1 
dt 

a 

A 

^a+1 ■ 

t . 

a- 1 


a = 2 . . . • M-1 . 


The linear interpolation to the point sets can be symbolically written 
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where Lin denotes a linear interpolation. 

The outside (top) boundary is a circular arc defined by 

x° = -R cos 0 + X 

0 

/NQ 

y = R sin 0 + y^. 

The physical domain in two dimensions is the region between the 
inside and outside boundaries. A rectangular region can be mapped onto 
the physical domain based on information from the boundaries and a 
linear connecting function. The transformation is 

X = x^n + x^(l - n) 


T 

y = y n + y^{l - n). 
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j=m 

Given the computational grid 

J 1 
k=l 


the physical grid is 


j=n 

‘"Jk’ 

k=l 


where 


Xjk ' *k"j 




i'jk = * ^k<’ - 


6, = (0 - 6, . )c, + 8 . , 

k ' max min'^k mm 


tk = t„?k 0 1 ? 1 1 


Concentration of the grid near the inside boundary and in the nose- 
shoulder region is accomplished by intermediate transformations 


e'^n - 1 


e'- 1 


, 0 £ n £ 1 , 


and 


. 1 

(e^l - 1)0 + Be - Bc^) 


0 < ^ < 1 . 
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The parameters k, , and B govern the concentration of the grid 
points. The above derivation for the transformation between a computa- 
tion domain and a physical domain is for a two-dimensional slice of a 
flow field. A three-dimensional representation is 

N-jk " yjk<='"f5i<*max - *niin> ^ 

^ijk " ^^max “ ^min^ * 

^ijk " ^jk 
0 < 5 1 1 • 

The elements of the Jacobian matrix are: 

M = (I _ A . )y(cos[£(^ - ^ ^ 

8? ^^max ^min'"^' ‘■^'^max ^min' ^min-^^’ 

= [sin(C[^ - ^ . ] + ^ • )] , 

8n ’■ '^‘■^max ^min-* ^tnin'-^ 3n 

P- = [sin(?[^ - ? • ] + ^ • )] P , 

3? ^^‘-^max ^min-^ ^min'-^ 3<; ’ 


( 2 . 18 a) 

( 2 . 18 b) 

( 2 . 18 c) 
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^ = 
3n 




(|) • ] + 
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- CcOSfCCd) ”<!)•]■*■ 

'- '^‘-^max ^mirv^ niin 
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3z 
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/vT 


3x - ax'" 96 X /I TTx 9x^ at a? 

^ ^ IT ^ ^ ^ ^ ^ ’ 
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^ = n ^ ^ ^ + (1 - n) ^ ^ 
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Table 3. Data description for a one and one-half inch spike-nosed body 


Inside boundary 

— - -- -- — - ■ - ■ 1 


Pt 

X 

y 

Pt 

X 

y 

Pt 

X 

y 

1 

0.0 

0.0 

21 

.01042 

.01804 


41 

.12475 

.02357 

2 

.00003 

.00109 

22 

.01138 

.01856 


42 

.12494 

.02428 

3 

.00011 

.00218 

23 

.01236 

.01903 


43 

.12500 

.02500 

4 

.00026 

.00326 

24 

.01337 

.01945 


44 

.12500 

.02667 

5 

.00046 

.00433 

25 

.01440 

.01981 


45 

.12500 

.07417 

6 

.00071 

.00539 

26 

.01544 

.02012 


46 

.12500 

.07533 

7 

.00102 

.00644 

27 

.01650 

.02038 


47 

.12506 

.07606 

8 

.00138 

.00747 

28 

.01757 

.02058 


48 

.12525 

.07676 

9 

.00180 

.00847 

29 

.01866 

. 02072 


49 

.12556 

.07742 

10 

.00227 

.00946 

30 

.01974 

.02080 


50 

.12597 

.07801 

11 

.00279 

.01042 

31 

.02083 

.02083 


51 

.12649 

.07853 

12 

.00336 

.01135 

32 

.02163 

. 02083 


52 

.12708 

.07894 

13 

.00398 

.01225 

33 

.12000 

.02083 


53 

.12774 

.07925 

14 

.00464 1 

.01311 

34 

.12083 

.02083 


54 

.12844 

.07944 

15 

.00535 1 

.01394 

35 

.12156 

.02090 


55 

.12873 

.07948 

16 

.00610 ; 

.01473 

36 

.12226 

.02108 


56 

.12883 

.07949 

17 

. 00689 

.01548 

37 

.12292 

.02139 


57 

.25000 

.09929 

18 

.00772 

.01619 

38 

.12351 

.02181 

Out 

Side Boundary 

1 


19 

. 00859 

.01685 

39 

.12403 

.02232 


= .43333 ft. R 

= .52208 ft. 

Of = 69.5° 







k2 

= .0001 yo 

= -.07500 ft. 

6 q = 8.295° 


. 00949 

.01747 

; 40 

.12444 

.02292 

kl 

= 2.2 B 

= 2 
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Table 2. Data description for a one-half inch spike-nosed body 



Inside Boundary 

Pt 

X 

y 

Pt 

X 

y 

Pt 

X 

y 

1 

0.0 

0.0 

21 

.01042 

.01804 

41 

.04142 

.02397 

2 

.00003 

.00109 

22 

.01138 

.01856 

42 

.04161 

. 02428 

3 

.00011 

.00218 

23 

.02136 

.01903 

43 

.04167 

.02500 

4 

.00026 

.00326 

24 

.01 337 

.01945 

44 

.04167 

.02667 

5 

.00046 

.00433 

25 

.01440 

.07981 

45 

.04167 

.07417 

6 

. 00071 

.00539 

26 

.01544 

.02012 

46 

.04167 

.07533 

7 

.00102 

.00644 

27 

.01650 

.02038 

47 

.04173 

.07606 

8 

.00138 

.00747 

28 

.01757 

.02058 

48 

.04192 

.07676 

9 

. 001 80 

. 00847 

29 

.01866 

.02072 

49 

.04223 

.07742 

10 

.00227 

. 00946 

30 

.01974 

.02080 

50 

.04264 

.07801 

11 

.00279 

.01042 

31 

.02083 

.02083 

51 

.04316 

.07853 

12 

.00336 

.01135 

32 

.02163 

. 02083 

52 

.04375 

.07894 

13 

.00398 

.01225 

33 

.03666 

. 02083 

53 

.044^1 

.07925 

14 

.00464 

.01311 

34 

.03749 

.02083 

54 

.04511 

. 07944 

15 

.00535 

.01394 

35 

.03822 

.02090 

55 

.04550 

.07948 

16 

.00610 

.01473 

36 

. 03892 

.02108 

56 

.04549 

.07949 

17 

.00689 

.01548 

37 

.03958 

.02139 

57 

.16687 

.09929 

18 

.00772 

.01619 

38 

.04018 

.02181 





19 

.00859 

.01685 

39 

.04070 

.02232 

Outside Boundary 

Xq = .4333 ft. Radius = .52208 ft. 

Of = 51.55° 

20 

. 00949 

.01147 

40 

.04111 

.02292 

\<2 

-0001 Vo = .07500 ft. 

2.2 B - L 

^^0 ^ 8.295° 
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Grids generated with this application of the "two-boundary technique" 
to the bodies shown in Figures 9 and 10 are shown in Figures 11 and 12. 
Tables 2 and 3 give the data used to generate the grids. The use of 
the grid in the Navier-Stokes solver is discussed in Chapter 4. 

2.5 Initial and Boundary Conditions 
Initial conditions are free stream conditions except at solid 
boundaries where no slip is imposed on the velocity. Free stream con- 
ditions are established from the Mach number = M^o, Reynolds num- 
ber E R = Re, characteristic length e L, and free stream tempera- 

“oo 

ture E T^. The speed of sound is 

(y(y - 



G 2 


■ft* 2 1 

where y = 1.4 and C„ = 4290 for air which is considered 

V 'sec- deg 

a perfect gas. The free stream velocity is 


and 


Uco = Moo Coo , 


= 0 . 


w = 0 . 

CO 


The free stream viscosity is 


2.27 X lO'^T 


3/2 


T + 198.6 

00 


The free stream density, pressure, and energy are 


Poo = Pco^^ooL/Re > 

oo 

P = pRT , 

00 oo ' 

and 


= CyT„ 


OO 
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No slip boundary conditions are imposed on the velocity at solid 
walls. A solid wall is considered to be isothermal and the tempera- 
ture Ty is fixed. The boundary condition on energy is e^^ = C^Tj^. 

The solid wall boundary condition for density = is obtained 
through a condition on pressure at the wall and the relation of density 
to pressure in the equation of state. The wall pressure boundary con- 
dition is obtained by approximately satisfying the momentum equations 

at the wall. Assuming that the gradient of the shear stress is zero 

9P 

at a solid wall implies that "^=0 where N indicates the normal 
direction. 

In general, the zero pressure gradient boundary condition at a 
solid surface can be enforced given the direction cosines Yy» Y2) 

of the normal vector on the surface. Then 


9N 



la 

9X 

la 

9y 

la 

3Z 



(2.19) 
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2.5.1 Boundary Conditions for Supersonic Flow About 
Wedqe-Cyl i nder Corners 

For wedge-cylinder geometries the transformation between the com- 
putational domain and the physical domain is given by Equation (2.17). 
The upstream boundary conditions at ^ = 0 are the free stream condi- 
tions. Solid walls occur at n = 0 and ? = 0 for c > 0. For n = 0 

the condition = 0 -»■ ^ = 0 where r is the radial direction from 
9N 3r 

the center line of the cylinder. 


1P| -lE.iZ+lE.9z + lPM = 

3r|^_Q ~ 3y 3r 3z 3r 3x 3r 

From Equation (2.17) 

1^ = cos (?02 + (1 - ^0]) = cos e 

1^ = sin (^02 + (1 - ^0^) = sin 0 

— = 0 
9r 

Therefore, 


3P 

3r 


n=0 


f cos e.f Sine 


and 


3r 


H=0 


,3P 3C . 3P 3n . 3P S 

%3y 3y' ® 
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where 


ac 

3? 

9y 

3z 

3n 

3n 

3y 

3Z 


3C 

9y 

3Z 


E elements of the Jacobian matrix at n = 0. 


Since gy 22 ~ 0* 


3P 

3r 


n=0 


(■1^ 1^ + 1^) cos 0 + (|^ + 1^ 1^) sin 6 

'3n 8y 3? 3y^ '3n 9Z 3? 3z^ 


= 0 . 


3P . 


3n 


is approximated by the one-sided difference 


3P ^ ^^l.k ~ ^^2,k ^3,k 

2An 


9P 

and, — approximated by the central difference 

oQ 


IP , '*2.l;-H ~ ’’z.k-l 
■ 2A5 


Then 


3r 


^^1 u “ ^^2 k ^3 k ^2 k +1 ~ ^2 k 1 

= ( — — — C, + r = 0, 

n=0 2An ' 2AC 


where 


Cl = ly cos 6 + Ij sin 6 


^2 " ly ® If ®- 
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Consequently, the pressure on the boundary n = 0 is approximated by 


l,k 


^^2.k ' ^3,k 


— fP - P ) 

3A? (P2,k+1 P2,k-1^- 


The normal pressure gradient boundary condition 
the wedge surface of a wedge cylinder corner is 


8P 


0 for 


C=0 


Iff ' H ♦ H> + Iff (|J '*’> * Iz 

■'ll <-='■"♦>■' If - *'> = ° 

where the directional cosines are: “ sin (p, y^ = 0, y^ = cos (p. 

The finite difference approximation for P at ? = 0 is 

^^J.2 ~ ^,j,3 _ ^ 


^,j,l 


where, 


3An D 


(P. ..1 9 - P 


1 


i ,j+l ,2 ■ i ,j-l ,2 


) 


2A? D2 

3A5 




D, = II (-sin ♦) + ||cos 4.. 
'’2 ' lx Iz 


D3 = -s1n 4. If . 


and 
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The far field boundary conditions (? = 1, n = 1» C = 1) are 


9U 

9? 


^=1 


= 0 ^ 
3n 


n=l 


= 0 . 


and 


?=1 


= 0. The condition fl" = 0 


implies that there is no change in the state variables with a change 
in C* The approximation is " ^i-T condition = 0 implies 
that there is no change in the state variables with a change in y. For 
the wedge-cylinder corners, this implies two-dimensional flow on a flat 
plate and/or inclined plate. That is: 


^=iy.^+M9n + 9Uii = 

9y 9? 9y 9n 9y 9c 9y 


Noting ^ ~ 0 and applying the finite difference approximation 


M _ ,^ 3 ,k ~ ^.1-1, k^ 9n ^ /^--l,k+l ■ ^j-l,k-lx 9c _ n 

3y ^ ^ Sv/ ^ ' / '577 “ ^ 


An 


9y 


2AC 


9y 


implies 


I = U + 

^j,k ^--l,k^2AC 


K 

fii 

M ^ j-1,k+l 
9y 


Vi,k-1> 


The condition ^=0 implies that there is no change in the state 
variables at c = 1 or that symmetry is imposed. In either case 
Ck = 1 for this study. 
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2.5.2 Boundary Conditions for Supersonic Flow 
About Spike-Nosed Bodies 

The transformation between a computational domain and the physical 
domain for spike-nosed bodies is given by Equation (2.18). In this 
case the solid surface is at n = 0, 0 < ? < 1, and 0 < ^ < 1. No 
slip conditions are imposed on the velocity and the temperature is 
fixed. The pressure is found in a similar manner to that used for the 
wedge-cylinder geometries except a less accurate approximation to the 
zero pressure gradient is used. In this case a first order approxima- 
tion is used and the pressure on the solid boundary is set equal to the 
pressure at the grid point next to the boundary. That is: 

^M,k " ^,2.k- 

The density at the boundary is found by the application of the equation 
of state using the boundary temperature and computed pressure. 

At c = 0, 0 < n < 1 , and 0 < C < 1 a symmetry boundary condi- 
tion is imposed. The three-dimensional grid is obtained by rotating 
a two-dimensional slice of the grid about the axis of symmetry. The 
line c=0, 0<n<l, and 0 < ? < 1. is coincident with the line of 

symmetry. The Jacobian matrix (Eq. (2.7)) is singular along this line. 
This does not create a problem, however, since the condition is imposed 
without using the transformation data from this line. The symmetry 
condition is 
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3U 


C=0 


= 0 or U 


i.j.l 


= u.. 




Free stream boundary conditions are imposed at n = 1 » 0 < c < 1, 

and 0<S<1. At <;=1. 0<ri<l»and 0<5<! a no-change 
boundary condition is imposed. That is: 




C=1 




3. COMPUTATIONAL ASPECTS 

The computational requirements for computing a grid and the 
Jacobian matrix associated with a grid using the "two-boundary tech- 
nique" are relatively minimal. It is however, necessary to plot the 
grid to visually assure that the desired constraints are satisfied. 
Ultimately, this phase of problem solving should be in an interactive 
mode with high bandwidth communications between the computer and a 
graphics terminal. The transformation data for an acceptable grid can 
be stored on a permanent file for later use. An alternate approach is 
to program the equations for the Jacobian matrix within a program for 
the solution of a flow-field. In this manner only parameters for 
the grid generation must be supplied. 

The computational requirements for the solution of the three- 
dimensional Navier-Stokes equations are extreme and tax the capability 
of any presently existing computer [1]. The approach taken in this 
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study is to adapt a viable numerical technique to the available large 
scale computer. The computer is the STAR-100 or its successor the 
CYBER 203. The computer architecture is based on vector processing 
with virtual memory storage. In addition to the vector architecture 
there are two aspects relative to the computer that have been very 
important in this study; (1) the capability of halfword arithmetic 
and storage; and (2) the effect of data transfer between secondary 
memory and primary memory. Halfword arithmetic has been used almost 
exclusively in the computations discussed later allowing for much 
larger grids than would be otherwise possible. Frequent transfers of 
data from primary memory to secondary memory and back have been 
minimized or avoided by constraining the grid size on which a solution 
is attempted. The transfer of data to and from secondary memory is 
relatively inefficient and is discouraged by high cost to the user. 

Another computational aspect is that the "Navier-Stokes solver" 
is relatively general. The application of initial and boundary condi- 
tions are performed in separate subroutines from the general solution 
procedure. Defining a new problem by initial and boundary condition 
does not require major programing modifications. 

In this chapter the MacCormack time-split algorithm is examined 
and vectorized for the CYBER 203 computer. Also, the program organiza- 
tion and how it relates to the virtual memory is presented. 
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3.1 Computational Technique 

The computational technique used in this investigation is the 
MacCormack time-split predictor-corrector algorithm [2-4] which was 
proposed about 1970 and is a derivative of the MacCormack unsplit pre- 
dictor corrector algorithm [26]. Both techniques are explicit which 
implies that they are time step stability limited [27]. Also, both 
techniques are second order accurate, and many investigators have been 
highly successful in applying them to a variety of fluid flow simula- 
tions [1,3,4,28]. An advantage of the MacCormack techniques is that 
they are relatively easy to apply to the transformed equations of 
motion (Eq. (2.3)). The split operator algorithm has the added 
advantage that different time step magnitudes can be used in each 
operator. A third hybrid scheme [29] can be applied by subdividing 
the operators into implicit and explicit portions. This approach, how- 
ever, is more complex and the success of its use is somewhat case 
dependent [30-31]. The explicit time-split technique has been chosen 
for this investigation because of its simplicity and vectorization 
characteristics for application on the CYBER 203 computer, however, 
both the unsplit and split algorithms are presented herein for contrast 
and clarity. 

3.1.1 MacCormack Technique 

The unsplit algorithm has the following two steps applied to the 
transformed equations of motion (Eq. (2.3)). 
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Predictor step : 


i > j » k 
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n 

i.j.k 
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j.k 


An 


'fj - fj-i) U J ^ <s- - 


^0-1 ^ 3y 
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Corrector step ; 


ljn+1 = i 


A? 




»i>i* 


j»k 


At 

An 






'i,k 


At 

AC 


(•'W - fk) If ^ " <Vl - \> t ^ + (H 


3C 


3y " '“k+1 


H ) — k 
”k^ 3z 


1 .J 
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This algorithm is applied for a time step by passing through a 
data base consisting of the state variables (p, pu, pv, pw, pe) and 
the transformation data and applying the predictor step. The corrector 
step is applied with the output of the predictor step, and the old 
state variables in the data base are replaced with the new values. 

The algorithm is repeated until a steady state solution is reached or 
an otherwise chosen stopping point is reached. 

3.1.2 MacCormack Time-Split Technique 

The split algorithm consists of a predictor and corrector step 
for each coordinate direction. Consequently, a predictor and corrector 
step for a coordinate direction is called an operator for that direc- 
tion (i.e., L.. (time step)). A time step is completed in 

direction ^ 

this algorithm with the application of each operator applied sym- 
metrically about the operator for the coordinate direction of primary 
flow. That is, for the corner flows studied herein 



At = At^ = ^ At^. 

n c 2 5 


where 
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For the spike-nosed body flow 


u"*l 1, 

1 .J.k 


K'^>] l'n<"‘c>l 




where 


At- = At = y At_. 

K n 2 c 


Each operator is defined by an output state solution for a 

given input state solution Therefore, 


W - 


where 


Predictor step : 
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in 

i J,k 





+ (Hi - Hi.-,) 
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Corrector step ; 


..out 


Atj 

j,k ^ AS 


<Fhi - i’ 
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-G.)|f 1 MH, -H..,)|fi] ] 

1 I/ / 


= C.k 


v/here 


Predictor step : 


U,. .. 
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At 


i.j.k i,j,k An 


<Fj - Fj-l) t j " 


'j "j-V 9y 


i.k 



Corrector step : 
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out 
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i-^(At^) = U 


out 
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Predictor step : 

i >j 
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Corrector step; 


I, out _ 1 /|,in . n- ^ 

2 “i.j.k " a? 


fc - F ) — k 


^^k+1 ■ ^k^ 3y ^^k+1 ■ ^k^ 3z ^ 


H 


i.J 


The unsplit algorithm requires only one pass through the data base per 
time step while the split algorithm requires several passes through the 
data base per time step. When the data base exceeds the primary stor- 
age capacity of the CYBER 203 computer a time penalty is imposed when 
data is called from secondary memory. However, a data management pro- 
cedure has been implemented to minimize the penalties associated with 
the use of secondary memory. 

It has been noted that the MacCormack algorithms are second order 
accurate. Forward and backward differences are applied such that after 
the predictor and corrector steps are completed an effective central 
difference approximation is obtained [3], This is demonstrated with 
derivatives of velocity components required in the viscous stress 
terms. Consider 


3u 

9 ?’ 


3n 


and 


3u 

9? 
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For a predictor step: 
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The unsplit and split MacCormack algorithms are time step stabil- 
ity limited, and there is no complete stability analysis to indicate 
the maximum allowable time step. A conservative time step employed by 
Shang and Hankey [3-4] has been used. This time step is 


At < min 


M + M + M 

Ax Ay Az 




-1 



where 


c E local speed of sound, 


A point that must be considered in the application of the 
MacCormack algorithms to viscous compressible flow with strong shock 
waves is the inclusion of terms to dampen oscillations in the region of 
a shock. A pressure dampening term suggested by MacCormack [31] is 
included in the finite difference approximation to the equations of 
motion. This term is 




36 , 


+ c 3 P 


4P 


36 


JlJ 


3U 

36„ 


Z = 1 ,2,3 


where 


<S-i = 


62 = n, and 6^ = ^. 
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3.2 Application of Vector Processing to the Computational Technique 
The MacCormack time-split algorithm has been programmed to run on 
the CDC STAR-100 and CYBER 203 computers. The program called the 
Navier-Stokes solver was first written in STAR FORTRAN and is described 
by Smith and Pitts [32]. The program has since been written in the 
SL/1 language [20] where 32-bit arithmetic is used to increase the com- 
putational speed and incore storage. 

3.2.1 Vector Processing Using the CYBER 203 Computer 

The CYBER 203 is a vector processing computer capable of achieving 
high result rates when a high degree of parallelism is present in the 
computation. When an identical operation is to be performed on consecu- 
tive elements in memory, a vector instruction is issued to perform the 
operation. Each vector instruction involves a time penalty, called 
vector startup, regardless of the length of the vector. As the length 
of the vector increases, the operation becomes more efficient since the 
penalty becomes relatively less important. 

The CYBER 203 has about one million words of primary memory with 
virtual memory architecture. Memory is referred to as pages. The two 
page sizes on the CYBER 203 are "small" pages which are 512 64-bit words 
and "large" pages which are 65536 words or 128 small pages. A user can 
have access to about 15 large pages in primary memory at any one time. 
The movement of data from secondary memory into primary memory involves 
moving pages of data in and out of primary memory. This is called a 
"page fault" and involves a startup time and transmission time just as 
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vector operations do. It then becomes important to make the most 
efficient use of the data when it is in primary memory in order to 
avoid a situation where the machine is spending more time moving pages 
of data in and out of primary memory than it is spending on actual 
computations. This is often referred to as "thrashing." Storing the 
data for a large data base program, such as a three-dimensional 
Navier-Stokes solver, in a conventional manner could very possibly 
lead to this situation. If, however, you design an interleaved 
data base [33] where the variables that are currently being used are 
stored together then it could result in less movement of "pages" of 
data. 

A capability of the CYBER 203 architecture is half-word arith- 
metic. This means that computation can be performed with 32-bit 
words and approximately two million 32-bit words can be stored in pri- 
mary memory. The speed of computation is approximately twice that 
achieved with 64-bit operations. A high level programming language 
SL/1 [20] is used to access the half-word capability. It is shown in 
the next chapter that for the MacCormack time-split algorithm 32-bit 
arithmetic is quite adequate and the computational speed and primary 
storage are approximately doubled as compared to FORTRAN version which 
uses 64-bit words. 

3.2.2 Program Organization and Data Management 

The Navier-Stokes solver is a derivative of a serial FORTRAN 
code which operates on the CDC-CYBER 175 and 7600 computers [2-3]. 
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Seven of the nine elements of the Jacobian matrix are programmed into 
the serial code and some redundant computation is performed to coexist 
with available memory. The Navier-Stokes solver is written in SL/1 and 
maintains the MAIN program logic and time step calculations found in 
the serial code (Fig. 13). The operator calculations are redesigned 
around the vector architecture of the CYBER 203 computer. Also, all 
nine elements of the Jacobian matrix are included in the new code. 

Normally it may be thought that vector lengths should be equal to 
the total number of grid points and vector operations sweep through the 
entire grid with each variable dimensioned to the number of grid points. 
However, for the number of variables involved and the large number of 
grid points this leads very quickly to the "thrashing" situation 
described earlier.. Instead, vectors are computed in planes in the 
^ directions (Fig. 14) with vector lengths approximately equal to the 
number of grid points in a plane. Temporary reusable vectors are main- 
tained for three local planes and a four-dimensional array S(I,L,J,K) 
contains the five state variables and nine elements of the metric 
coefficient for each grid point. 

Forming S(I,L,J,K) is the essence of interleaving the data base. 
The index L refers to the 14 variables. The index I refers to the 
plane in the ^ direction and the indices J and K refer to the 
grid points in a plane. In this manner all variables for each plane 
are stored in contiguous locations. Computation proceeds from the 
first plane to the last plane for each operator. 
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In order to minimize the sweeps through the data base S(I,L,J,K), 
the corrector step for a plane is performed as soon as enough planes of 
the predictor step are available (Fig. 15). Consequently for the 
application of each operator, there is one sweep through the data base 
and five sweeps for a time step. 

Within a C plane, a vector sweep is from the lower left hand 
corner to the upper right hand corner. The exact length, starting 
point, and end point of the vector is dependent on the operator and the 
direction of the finite differencing. Vector operations include 
boundary points where erroneous values are computed during a vector 
computation. The boundary condition subroutine is called to compute the 
boundary conditions and overwrite the erroneous values. 

The transformation data which consist of the nine derivatives of 
the computational coordinates with respect to the physical coordinates 
at each grid point are computed in a separate program and stored on a 
disk file. Once a geometry is established, the transformation data 
remains constant. It is read from the disk file by the Navier-Stokes 
solver and stored in the S array (L = 6... 14). 

Externally the Navier-Stokes solver operates like the serial pro- 
gram starting with some initial conditions and integrating with the 
finite-difference algorithm until "steady state" is reached, or the 
program is stopped. A restart capability is included so that a solu- 
tion can be obtained in several runs with intermediate observation of 


the solution. 
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4. RESULTS AND DISCUSSION 

Assuming that the equations of'motion (Eq. (2.3)) are valid and 
that a finite difference technique (MacCormack time-split algorithm) is 
used to numerically solve them, the physical grid and its relation to 
the computational grid form the foundation for a solution. This chapter 
is aimed at establishing the applicability of the "two-boundary tech- 
nique" for grid generation by obtaining solutions to complex flow 
fields using transformation data derived from the technique. Also, the 
robustness of the MacCormack technique and the computational capabili- 
ties of the CYBER 203 computer are demonstrated. 

The "two-boundary technique" is applied to two distinctly differ- 
ent supersonic flow problems. In Chapter 2, relationships between the 
computational grid and the physical grids for a family of three- 
dimensional corners and spike-nosed bodies are derived. The derivatives 
for the transformation data are also presented. In this chapter the 
transformation data and boundary conditions for the flows are applied 
in the Navier-Stokes solver. The three-dimensional corner flow 
fields are extremely complex with shock-boundary layer interactions and 
three-dimensional separation. The grids are concentrated near the solid 
boundaries and in the corner. The flow about the spike-nosed bodies is 
characterized by a strong bow shock and a highly separated region 
between the nose and shoulder. The grid is concentrated in this region, 
and the solutions are unsteady. The two flow problems are considered 
separately, and the grid characteristics are pointed out. 
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4.1 Supersonic Corner Flow Using Two-Boundary Grid Generation 

Supersonic flow about three-dimensional corners occurs in many 
high speed aerodynamic situations. Over the years there have been many 
experiments to study this flow phenomenon [34-38]. More recently there 
have been numerical experiments to compute supersonic flow about three- 
dimensional corners. Inviscid compressible solutions about planar 
three-dimensional corners have been obtained by Kutler, Shanker 
et al . [39-40] and Marconi [41]. Inviscid solutions, however, do not 
account for the strong inviscid-viscid interactions which occur. 
Asymptotic viscous solutions have been obtained by Weinberg and 
Rubin [42] and Ghia and Davis [43]. These solutions require extensive 
assumptions about the flow and also do not adequately describe the 
inviscid-viscid interactions. 

The solution of the compressible Navier-Stokes equations is the 
most conclusive way to compute supersonic flow about three-dimensional 
corners, and there are several published solutions. Shang and Hankey [3] 
using a time split MacCormack technique for the solution of the Navier- 
Stokes equations simulated the flow about an asymmetric 15° wedge-plate 
corner at Mach number 12.5. This numerical simulation corresponds to 
the physical experiment of Cooper and Hankey [38]. For this simulation 
the largest computational mesh was 8 x 32 x 36. In a later numerical 
study Shang and Hankey [4] computed a turbulent flow about a sym- 
metric 9.48° planar corner at Mach number 3. 

Hung and MacCormack [44] computed supersonic laminar flow about 
a 10° asymmetric planar corner preceded by a rectangular corner. They 
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used a time-split MacCormack technique where the conventional two-step 
scheme is used in the operator and in the Ly and operators 
in the regions far from the solid boundaries. In the inner regions 
the Ly and operators are further split into hyperbolic and 
parabolic operators for the inviscid and viscous terms. Later Hung 
and MacCormack [45] extended the code for the accelerated technique to 
compute turbulent flow, and Horstman and Hung [46] studied several 
planar turbulent corner flows. 

Most published discussion on supersonic flow about three-dimensional 
corners including that mentioned above has dealt with planar intersecting 
boundaries. This means that rectangular coordinates or a simple trans- 
formation have been used. Herein, both planar intersecting and planar- 
cylinder intersecting corners are discussed. In Chapter 2, Equa- 
tion (2.17), derived using the "two-boundary technique," relates a 
physical domain with wedge-cylinder boundaries to a computational 
domain. A plane-cylinder corner is formed by letting the wedge angle 
be zero and planar intersecting corners are approximated by letting the 
radii be very large. The transformation data obtained by differen- 
tiating Equation (2.17) is used in the Navier-Stokes solver along with 
the boundary conditions discussed in Chapter 2. 

During this investigation many solutions with varying Mach number, 
Reynolds number, and grid size have been obtained. Portions of the 
following discussion have been presented by the author in [47]. The 
initial numerical experiment is the solution of the flow about a rec- 
tangular corner (Fig. 5a) with a 31 x 31 x 31 grid. The velocity 
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solution u/u„o at x/xl = .554 where the flow is two-dimensional is 
shown in Figure 16 and is compared with that presented by Hung and 
MaCCormack [44j. The Mach number is 3 and the Reynolds number is 
2.78 X 10^/m (7 x 10^/in.). The concentration of the grid based on 
Equation (2.17a) has a value of k-j = k£ = 3.8. The agreement of the 
velocity with that of the reference is good. This initial solution is 
obtained on the STAR- 100 computer with the FORTRAN version of the 
Navier-Stokes solver. The computational rate is 1.5 x 10"^ seconds per 
grid point per time step. The remaining corner flow solutions to be 
discussed have been obtained on the CYBER 203 computer and the SL/1 
version of the Navier-Stokes solver. The computational rate for the 
largest grid used is 4 x 10"^ seconds per grid point per time step. No 
significant degradation of the solutions using the small 32-bit word 
size observed. 

The next step in this experiment is to obtain the solution of a 
family of wedge- cylinder and wedge-plate corner flows using a 
20 X 36 X 36 grid. There are three planar corners with 0, 6, and 
12.2° wedge angles and three wedge- cylinder corners with the same 
wedge angles. Equation (2.17) has been used to define the grids and 
transformation data. The physical dimensions of the domain are shown 
in Figure 17, and the grid at x/X|^ = 1 are shown in Figure 7. The 
concentration parameter in Equation (2.17a) is 3.8. Two additional 
solutions for a planar 18° corner have been obtained — one with the 
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Fig. 17 Physical dimensions of corners. 
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Fig. 18 Three-dimensional corner flow 
characteristics. 
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Reynolds number used in all other solutions, and the other with a 
Reynolds number equal to 3.9 x 10^/tfi (1 x 10^/in.). The physical posi- 
tion of the right boundary is moved further to the right because the 
effects of the shock are too close to the original position. The Mach 
number is 3.64 and the Reynolds number except for the one 18° planar 
wedge corner is 2.92 x lO^/m (7.42 x 10^/in.). The free stream and 
body temperatures are respectively 217 K (390 R) and 607 K (1092 R). 

The fluid properties are for air. Solutions are started impulsively 
except for the one high Reynolds number solution and are run for 
approximately eight characteristic times for a steady state solution. 
For the high Reynolds number case, the initial state is that of the 
lower Reynolds number solution. A steady state is obtained through 
several increments until the higher Reynolds number is reached. 

Flow visualization is presented in the form of continuous tone 
density distributions and a combination of continuous tones and con- 
tours for the temperature at x/x^ = 0.83. These distributions are 
made for this document from color distributions obtained on the Dicomed 
Graphics System mentioned in the introduction. Velocity vectors (uxw) 
at the first grid point above the flat plate and cylinder surface are 
shown to indicate the direction of flow near the base surface. 

Velocity vectors (vxw) at x/x^ = .83 are shown to indicate the cross - 
flow velocity field. 

The primary observation from this experiment and others is that 
the flow is basically conical. Also, the wedge-cylinder flow fields 
for the cylinder radius used are very similar to the planar corners 
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for the same wedge angle. For the 0 and 6° wedge angle solutions, the 
zones of flow (Fig. 18) described by Charwat and Redekopp [34] are 
readily seen (Figs. 19, 20, 21, 22, and 23). That is, for a traverse 
plane of a symmetric planar corner, there is Zone I (a region of 
conical flow bounded by slip lines and the corner shock). Zone II (a 
region of complex flow bounded by slip lines and a strong internal 
shock). Zone III (an outer interaction region characterized by a com- 
pressive fan centered at the triple shock intersection point), and 
Zone IV (undisturbed wedge flow). Figure 18 depicts the flow situation 
in a symmetric traverse plane. Charwat and Redekopp [34] also concluded 
that the flow structure remains qualitatively similar with change in 
Mach number for symmetric planar corners, and is distorted without 
losing its identity for asymmetric planar corners. The rectangular 
corner solution (Fig. 19) shows the zones of flow in a symmetric 
pattern. The plate-cylinder corner solution (Fig. 22) is very similar 
to the rectangular corner flow. The base surface curvature modifies 
the symmetry to some extent, but close to the corner the flow is 
almost identical. The temperature distributions for the rectangular 
corner and plate-cylinder are very orderly and follow the zones of flow. 

As predicted by Charwat and Redokopp [34] the 6° wedge corners 
show the zones of flow somewhat distorted (Figs. 20 and 23), Again the 
6° wedge-cylinder corner solution is very similar to the 6° planar 
corner solution. Also, the 6° corner solutions show two interesting 
characteristics that are observed in the larger wedge angle solutions 



but not observed for the 0° solutions. The crossflow velocity 
separates on the wedge surface about the internal shock. 

The crossflow velocity separation is observed in all the non-rzero 
wedge angle cases computed at Reynolds number equal to 2.92 x 10^/m 
(Figs. 10, 21, 23, 24, and 25). The second characteristic is the high 
crossflow velocity near the flat plate or cylinder surface under the 
internal shock. This phenomenon is observed on all the non-zero wedge 
angle solutions. 

Cooper and Hankey [38] performed experiments on a 15° asymmetric 
planar corner at Mach 12.5. They observed one triple point instead of 
the two observed by Charwat and Redekopp [34]. The numberical experi- 
ment by Shang and Hankey [3] basically confirmed the one triple point 
observation. The 12.2 and 18° wedge angle solutions obtained at the 
low Reynolds number and described herein appear to show two triple 
points although highly distorted. The corner shocks are almost 
vertical. The high Reynolds number solution for the 18° wedge angle, 
however, indicates that there is only one triple point. In fact, the 
high Reynolds number solution is qualitatively very similar to that of 
Shang and Hankey [3]. The relevance of Reynolds number is discussed 
further at a later point. 

Korkegi [48] has described supersonic flow in three-dimensional 
corners and associates separation with the shock strength which is a 
function of the wedge angle. Both laminar and turbulent corner flow 
display the separation phenomenon, however, a larger shock strength is 
required to cause separation in a turbulent flow. Korkegi [48] 
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describes three-dimensional separation as a line along which the flow 
lifts off a continuous solid surface and three-dimensional reattachment 
as a line of flow impingement on a continuous solid surface. Further, 
secondary separation can occur when the shock strength is sifficient 
for the reverse flow region to separate. The convergence and divergence 
of the uxw velocity vectors close to the surface indicate separation 
and reattachment. The 0 and 6° wedge angle corner solutions (Figs. 19, 
20, 21, and 23) show no strong indication of separation on the plate or 
cylinder surfaces. The uxw velocity vectors align themselves with 
the direction of the wedge inside the shock and align themselves some- 
what with shock in the shock region. The 12.2° wedge corner solution 
(Figs. 21 and 24) show evidence of convergence of the uxw velocity 
vectors or separation, but divergence is not clear. 

The 18° wedge corner solution at the lower Reynolds number 
(Fig. 25) displays both the convergence and divergence of the uxw 
velocity vectors. Figure 26 shows the solution of the 18° wedge corner 
at a Reynolds number of 3.9 x 10®/m (lx 10^/in.). The plate surface 
pressure (Fig. 27) shows the trough-like pressure variation described 
by Korkegi [48j as an indication of secondary separation. The 
interesting point about the 18° wedge corner is the effect of 
increasing the Reynolds number. Starting with the steady state solu- 
tion for the low Reynolds number 18° wedge corner solution, the 
Reynolds number is gradually increased and a new steady state achieved. 
The velocity vectors are rescaled since the velocity is much greater 
in magnitude at the same relative position as the low Reynolds number 
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Fig. 19a Density distribution Fig. 19b Temperature 


Re = 2.92 X 10^/m 



Fig. 19c u X w Velocity Fig. 19d v x w Crossflow velocity 


Fig. 19 Flow field description for a rectangular corner. 
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Fig. 20b Temperature 


Fig. 20a Density distribution 



Fig. 20 Flow field description for a 6° wedge-plate corner 




Fig. 22a Density distribution Fig. 22b Temperature 


Re = 2.92 X 10^/m 



Fig. 22c u x w Velocity Fig. 22d v x w Crossflow velocity 


Fig. 22 Flow field description for a plate-cylinder corner. 
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Fig. 23a Density distribution 


Fig. 23b Temperature 


Re = 2.92 X 10^/m 



Fig. 
Fig. 23 


23c u X w Velocity Fig. 23d v x w velocity 

Flow field description for a 6*^ wedge-cylinder corner. 
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Fig. 24a Density distribution 


Fig. 24b Temperature 
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Fig. 24c u x w Velocity Fig. 24d v x w Crossflow velocity 


Fig. 24 Flow field description for a 12.2 wedge-cylinder corner 
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case. The grid is coarse for defining the boundary layer, but the 
solution correlates better with physical experiments than the low 
Reynolds number solutions. Only one triple point is observed and 
separation, reattachment, and secondary separation are readily seen in 
Figure 26 and the surface pressure in Figure 27. Obviously the Reynolds 
number is a key parameter in the inviscid-viscid interaction in a three- 
dimensional corner. 

From a qualitative point of view it is apparent that the grids 
generated with the "two-boundary technique" perform well for the low 
Reynolds number solutions. The surface pressure for the 0, 6, and 12.2° 
wedge-plate and wedge-cylinder solutions are shown in Figures 28-29. 

The surface pressure for the 12.2° wedge-plate corner is later compared 
with a solution obtained for the same boundary geometry with a finer 
grid. 

4.1.1 High Resolution Grid Solutions 

It is implied in the above section that more grid points are 
needed to resolve the fine structure of supersonic flow about three- 
dimensional corners at high Reynolds numbers. In this section three 
solutions at a high Reynolds number and one at a low Reynolds number 
are described where a 12 x 64 x 64 grid is used. The objective of the 
section is to present the most highly defined solutions within existing 
computational capabilities and further validate the use of the "two- 
boundary technique." The coordinate transformation presented by Shang 
and Hankey [4] is used to compute the solution about a 12.2° symmetric 
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Fig. 29 Surface pressure for wedge-cylinder corners. 
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wedge-wedge corner (Fig. 30) at Mach 3.64 and Reynolds number 
3.9 X 10^/m (1 X 10^/in.). The solution is compared with the experi- 
mental data obtained by Charwat and Redekopp [34]. Using the Shang- 
Hankey transformation which can also be derived with the "two-boundary 
technique" and letting one of the wedge angles be zero, the solution of 
the flow about a 12.2° asymmetric wedge corner is obtained with the 
12 X 64 X 64 grid. A 12 x 64 x 64 grid is generated using Equation 
(2.17) and a solution is again obtained for the flow about the 12.2° 
asymmetric corner. The procedure for obtaining solutions at the high 
Reynolds number is to first start the solution at a low Reynolds number 
and increment to the larger value. The flow conditions are changed on 
the upstream plane and are integrated downstream. 

The transformation used in [4] is 

5 = x/X|^ (3.1a) 

n = r iln[l + (e*^ - 1)(J - tan 6^) |-] (3.1b) 

k L 

C = - iln[l + (e*^ - l)(f - tan 9) 4 --] . (3.1c) 

R ^ ‘-1 

This transformation can be derived with the "two-boundary technique" 
where 


x-,(5,c) = x(?) = 4 \ 




Fig. 3Q Hypothetical syiranetric 


corner flow from Ref. 37. 
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y-|(C.rJ = 5 X|_ tan 6 |^ 

= C X|_(tan 6 |^ + ^ x^(tan e^)(l - 5 ) 

X2(C,c) = x(0 = C Xl 
y2(C,0 = £ X^Ctan 6 ^^ + Y^) 


Z2(C.^) = z-|(C,c) 



0 1C 1 1 
0 < c < 1 

Using a linear connecting function and rearranging terms yields 


X = S X^ 

y = 5 XJtan 6 „ + 
z = c X|^{tan e^, + z,_(§l^)) 

0 < n < 1 


(3.2a) 

(3.2b) 

(3.2c) 
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Equation (3.1) is the inverse of Equation (3.2). This transformation 
applies only to corners formed from planar intersecting boundaries, and 
the transformation data is derived directly by differentiating C» n, 
and c with respect of x, y, and z. 

The solution obtained with a 12 x 64 x 64 grid and Equation (3.1) 
for a 12.2° symmetric wedge corner is shown in Figure 31-34. The 
density and temperature distributions show the zones of flow described 
by Charwat and Redekopp [34]. In this solution, the slip lines are 
very well defined. The line contour plot (Fig. 32) and the surface 
distribution of the density (Fig. 33) also shows the zones of flow and 
the rapid change in density associated with the slip lines. The slip 
line definition is not as apparent for the rectangular corner solution 
shown in the previous section. No crossflow separation is observed, 
which is the case with the rectangular corner. The position of the 
shock structure from [34] is superimposed on the density contour plot 
(Fig. 32) and the agreement is good. In [34] the shock structure is 
not perfectly symmetric, and only the upper half is used for comparison. 
The surface pressure in Figure 34 shows the comparison of the Navier- 
Stokes solution with Charwat 's experiment [34]. The plateau, dip, and 
overshoot described by Korkegi [48] occur at nearly the same position. 
There is, however, some disagreement in the pressure magnitude. The 
Navier-Stokes solution is on the low side. There is pressure dispersion 
in Charwat's data [34] but the pressure magnitude is still larger than 
that obtained in the Navier-Stokes solution. Two possible sources for 
this difference are (1) the Reynolds number may not be exactly the same 
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in the numerical and physical experiments, or (2) the pressure damping 
in the numerical technique may be retarding the pressure magnitude. 

The plateau region (Fig. 34) is associated with the relatively low uxw 
velocity vectors in Figure 31 and the separated region. The separation 
line is under the plateau near the undisturbed flow and reattachment 
occurs at the high pressure rise and where there is corresponding large 
velocity vectors. 

The dip in the pressure curve is associated with secondary separa- 
tion as described by Korkegi [48]. This secondary separation region is 
shown in the density contour plot (Fig. 32) and there are closed con- 
tour lines which indicate vortical motion as described by Watson [37]. 
However, the author has not observed vortical motion in the velocity 
data. 

It is evident that the 12 x 64 x 64 grid derived from Equation 
(3.1) allows for the definition of the fine structure of the flow. The 
solution obtained with the Navier-Stokes solver is in good comparison 
with the corresponding experiment. The possible concentration of more 
grid points in the secondary separation region may be desirable but 
would not work well with Equation (3.1) because of its exponential 
characteristics . 

Following the solution for a symmetric planar corner flow, super- 
sonic flow about an asymmetric corner is obtained by letting the wedge 
angle in Equation (3.1) on the bottom surface be zero. The wedge 
angle 0^^ remains 12.2° and the 12 x 64 x 64 grid for this corner is 
generated with Equation (3.1). The Navier-Stokes solver is applied 
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with the transformation data starting impulsively with a low Reynolds 
number and incrementing to 3.9 x lO^/rn. The solution is shown in Fig- 
ure 35-38, but is not compared to a physical experiment. The solution 
obtained with a grid generated using Equation (2.17) is compared to 
this solution. 

For the asymmetric corner and the 12 x 64 x 64 grid only one 
triple point is observed. However, the internal shock from the wedge 
generates a very similar flow pattern as in the symmetric corner near 
the bottom surface. The pattern is, however, closer to the wedge sur- 
face since there is no corner shock. The crossflow separation is 
observed at the high Reynolds number but is relative weak compared to 
the lower Reynolds number solutions. 

Equation (2.17) is used to generate a 12 x 64 x 64 grid with the 
previously used physical diiffensions and contraction parameters 
k-| = k2 = 2.9. The transformation data for the grid is applied in the 
Navier-Stokes solver in the same manner as before. The surface pressure 
for the solution at Reynolds number 2.92 x 10^/m is shown in Figure 39 
and compared with that obtained with the 20 x 36 x 36 grid. The agree- 
ment is good where the grid is dense and the solution at Reynolds 
number 3.9 x 10®/m is obtained through a series of increments. This 
solution is presented in Figures 40-43. The surface pressure in 
Figure 43 is compared to that obtained with the grid defined with 
Equation (3.1). Overall the agreement is good. One point that is 
noticed, however, is that the shock is more smeared in the solution 
obtained with the transformation data from Equation (2.17). This is 
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Fig. 37 Perspective view of distribution of density for 
12.2° asymmetric corner x/X[_ = .8 (12 x 64 x 64 grid - 
exact boundaries). 
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Fig. 38 Surface pressure for 12.2° asymmetric corner 
(12 X 64 X 64 grid - exact boundaries) . 
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Fig. 39 Surface pressure comparison for grid refinement 
at Re = 2.92 x loVni. 
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Fig. 42 Perspective view of distribution of density for 
12.2° assyinmetric corner x/xj_ = .8 (12 x 64 x 64 grid 
approximate boundaries). 
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attributed to the growth of the transverse mesh distribution in the 
windward direction. Using Equation (3.1) there is a very fine grid 
upstream and the grid becomes coarser with the expanding flow down- 
stream. Using Equation (2.17) the grid is coarse upstream and is finer 
at the downstream positions. 

It is concluded that the "two-boundary technique" in viable for 
computing grids about three-dimensional corners where the Navier-Stokes 
equations can be applied to compute supersonic flow. It is noted, how- 
ever, that if only planar intersecting corners are of interest. Equa- 
tion (2.17) is not the optimal application of the technique. Planar 
intersecting corners are approximated using Equation (2.17). Neverthe- 
less, the "two-boundary technique" is highly versatile and affords a 
great deal of flexibility. 

4.2 Supersonic Flow About Spike-Nosed Bodies 
Spike-nosed configuration occurs in many supersonic flow situa- 
tions. Unlike the three-dimensional corner flows that are studied 
herein, the flow fields about spike-nosed bodies can be highly unsteady. 
The unsteadiness manifests itself as self-sustained oscillations which 
have been observed for a wide variety of shear-layer impingement con- 
figurations [49]. Two spike-nosed bodies are examined in this section: 

(1) a body with a small nose length to should height ratio (0.71); and 

(2) a body with a large nose length to shoulder height ratio (2.14). 
These two configurations have been examined by Shang, Hankey, and 
Smith [50j. The first two authors posed the problem and this author 
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applied the "two-boundary technique" to generate the grids and obtained 
the initial solutions with the Navier-Stokes solver. The emphasis 
herein is on the application of the "two-boundary technique" for grid 
generation about the spike-nosed bodies. The "two-boundary technique" 
has been used to generate grids using a linear approximation to the 
body surface, a circular arc outer boundary, and a linear connecting 
function. A parabolic algebraic function and an exponential function 
are used to concentrate the grid in the nose shoulder region. The 
details of the application of the "two-boundary technique" to these 
geometries are developed in Chapter 2 along with the derivatives for 
the transformation data. The three-dimensional Navier-Stokes solver 
is used to obtain the numerical solutions to the flow fields. For 
these solutions the z-coordinate direction is the windward direction. 

An axisymmetric solution is obtained by the rotation of the grid about 
the z-axis and solving the three-dimensional equations of motion. The 
argument for this approach rather than developing an axisymmetric 
solver is the considerable time savings compared to the programming a 
specific code for this problem. These solutions demonstrate the versa- 
tility of the Navier-Stokes solver. 

4.2.1 One-Half-Inch Spike-Nosed Body 

The spike-nosed body with the small ratio of nose length to 
shoulder height has a nose length of 12.5 mm (0.5 in.) and a shoulder 
height of 19.05 mm (0.75 in.). The boundary surface for this body is 
shown in Figure 9 and the grid generated with the "two-boundary 
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technique" is shown in Figure 11. The data used to generate the grid 
is found in Table 2. Figure 44 shows the density distribution of the 
developing flow. A relatively low density region compared to the sur- 
roundings develops at the top of the nose and later sheds off the 
shoulder. After this point only low amplitude oscillations occur in 
the nose-shoulder region. Plots of pressure along three lines of the 
grid (Fig. 45) are shown after 4000 time steps in Figures 46-48. 

Two modifications to the data for generating the grid for the one- 
half-inch spike-nosed body are performed and the solutions recomputed. 
The objectives of the modifications are to assure the validity of the 
grid generation technique. 

The first modification is a change in the concentration of grid 
points. The new grid is obtained by making the constant k = 3 instead 
of 2.2 in Table 2. This means that there are fewer grid points to 
define the shock in the radial direction and more points near the inner 
boundary. Figures 49-51 show the pressure comparison with the original 
case at the 4000th time step for the three lines. The second modifica- 
tion moves the outer boundary closer to the inner boundary concentrating 
more points in the shock region. The radius of the circular arc defining 
the outer is changed from 13.26 mm (0.522 ft) to 12.55 mm (0.494 ft) and 
the constant k is set equal to 1.75 so that the distance of the 
nearest grid point to the inner boundary is approximately the same as 
that in the original case. Figures 52-54 show the density comparison 
with the original case and the first modification. It is seen from the 
plots that there is overall good agreement. The differences that exist 
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Grid lines for quantative comparison 
one-half inch spike-nosed body. 
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Fig. 47 Density along line twenty-nine for base case. 
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Fig. 49 Comparison of density solution for grid concentration 

change line = 1. 
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Fig. 50 Comparison of density solution for grid concentration 

change line = 29. 
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Fig. 51 Comparison of density solution for grid concentration 

change line = 53. 
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Fig. 52 Comparison of density solution for outer boundary change 1. 
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Fig. 54 Comparison of density solution for outer 
boundary change line = 53. 
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are due to the grid spacing in the shock region and inner boundary 
region. The shock is more smeared when there are fewer grid points in 
that region and the placing of grids points closer to the inner boundary 
affects the accuracy there. No attempt is made here to analyze the low 
amplitude oscillation for the one-half-inch spike-nosed body. It is 
seen that the perturbed grids produce good agreement with the original 
case and, this is an indication that the grid generation technique is 
viable for this application. 

4.2.2 One and One-Half-Inch Spike-Nosed Body 

The spike-nosed body with the large ratio of nose length to 
shoulder height has a nose length of 38.1 mm (1.5 in.) and a shoulder 
height of 19.05 mm (0.75 in.). The body surface for this configura- 
tion is shown in Figure 10 and the grid generated with the "two- 
boundary technique" is shown in Figure 12. Table 3 contains the data 
used to generate the grid. Table 3 is the same as Table 2 except the 
nose length is one inch longer and the outside circular arc is longer. 

The Navier-Stokes solution for this grid and initial conditions pro- 
duces high amplitude oscillation. The characteristics of the oscilla- 
tion can be seen from experimental observation. Figure 55 shows a 
sequence of shadowgraph pictures [51] of a spiked-nosed body flow 
field in one cycle of oscillation. After the initial transient, the 
bow shock interact with the shock from the shoulder near the top of the 
shoulder. A strong reverse flow occurs between the shock and the 
boundary forcing the shock to bulge out. With the shock bulged out. 
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the high energy fluid behind the shock along with shock is forced 
downstream by the outer flow until the bow shock again interacts with 
shoulder shock. The cycle repeats itself. This phenomena is also 
shown by Harney [52] for blunt nosed bodies which is approximated by 
the one and one-half-inch spike-nosed body. Widhopf [53] has shown 
the oscillating phenomena over indented nose tips while numerically 
solving the Navier-Stokes equations. Figure 56 shows a sequence of 
density distributions for one cycle taken from the numerical solution. 
The phenomenon described above is observed in the numerical solution of 
the flow about the one and one-half-inch spike-nose body. The cycle 
repeats itself and maintains its form. The oscillation phenomenon is 
independent of Reynolds number [50] but is dependent on the speed of 
sound which in turn depends on the free stream temperature. The pri- 
mary frequency is 3100 cycles/sec where the free temperature is 217 K 

fk 

(290 R). In [51] using a Reynolds number of 7.78 x 10 /m and free 
stream temperature of 111 K (200 R) a Tower primary frequency of 
2665 cycles/sec is observed. The pressure at approximately point 29 
(Fig. 48) and the observed pressure from an experiment conducted at the 
Wright-Patterson Flight Dynamics Laboratory are shown in Figure 57. It 
is seen in the figure that the wave form and frequency from the experi- 
ment and the numerical solution are in good agreement. 

The "two-boundary technique" for grid generation has been success- 
fully applied for two spike-nosed bodies. The Navier-Stokes solver has 
produced unsteady numerical solutions for several initial conditions 
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Figure 57.- Surface pressure on one and one-half inch spike-nosed body 
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with Reynolds numbers up to 7.78 x 10®/m. Overall, the solutions that 
have been obtained simulate the observed phenomena very well. 

5. CONCLUSIONS 

An algebraic grid generation technique has been developed and 
explored in conjunction with the solution of the compressible three- 
dimensional Navier-Stokes equations. The technique called the "two- 
boundary technique" is simple to understand, easy to apply, and has a 
high degree of generality for the finite difference solution of complex 
flow field problems. The "two-boundary technique" allows direct control 
of a grid and direct computation of the Jacobian derivatives. 

The viability of the grid generation technique is demonstrated 
.through the development and application of a Navier-Stokes solver which 
operates on the CDC CYBER-203 vector computer. The computer program is 
based on a MacCormack time-split technique which is chosen because of 
its compatibility with vector computer architecture. The finite differ- 
ence algorithm is written in the SL/1 programming language, and the 
32-bit word length arithmetic and storage option is used. This option 
doubles the number of grid points that can be used for a given amount 
of memory and approximately doubles the computational rate as compared 
to the normal 64-bit words. Using SL/1 and the halfword option the 
computational rate is 4 x 10“^ seconds per grid point per time step, 
and solutions with 5 x lo^ grid points can be obtained without using 
secondary memory. It is concluded from the numerical experiments 
presented in the present study that the 32-bit word length is adequate 
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when solving the Navier-Stokes equations for supersonic laminar flow 
using an explicit MacCormack technique. 

Complex supersonic flow field solutions are obtained for two dis- 
tinctly different geometries using the "two-boundary technique" for 
grid generation and the Navier-Stokes solver. First, supersonic flow 
solutions about a family of three-dimensional corners are obtained. 

These flow fields reach a steady state but are characterized by strong 
shocks and three-dimensional separation. The Mach number is 3.64, 
Reynolds numbers are 2.72 x 10 /m and 3.9 x 10 /m, and the fluid proper- 
ties are for air. It is shown that the solutions obtained agree well 
with physical experiments and other numerical experiments. Also, 
corner flow solutions with 5 x 10“^ grid points are among the most 
refined Navier-Stokes solutions obtained to date. The second flow 
situation is supersonic flow about spike-nosed bodies. In this case, 
the flow is axisymmetric, unsteady, and characterized by a strong bow 
shock and massive separation. The Mach number is 3 and the Renolds 
number is 7.78 x 10^/m. The numerical solutions show dramatically the 
oscillating flow generated by the interaction of the bow shock and 
shoulder wall of the body. The surface pressure and oscillation fre- 
quency compare very well with corresponding wind tunnel experiments. 

The successful numerical solution of the flow fields support the primary 
conclusion that the "two-boundary technique" is viable for generating 
grids for complex flow field solutions. Also, for the spike-nosed 
bodies, considerable development time for a specialized axisymmetric 
code is saved. 
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Plans for the use of the "two-boundary technique" include develop- 
ment of grids with wing-fuselage boundaries, analysis of non-orthogonal 
grids, development of additional spike-nosed body grids, and the 
development of numerical grid control functions. 
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